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Abstract 


A  paxticle-in-cell  ansatz  for  solving  the  Euler  equations  for  geophysical  fluid  dynam¬ 
ics  is  described.  The  approach  is  ideally  suited  for  “layered”  models  in  which  density 
and  velocity  are  independent  of  the  vertical  coordinate  in  fluid  layers  but  generally 
vary  with  layer  specification.  The  material  acceleration  terms  in  the  Euler  equations 
are  solved  at  each  particle  while  the  gradient  terms  are  evaluated  on  a  grid  and  inter¬ 
polated  at  each  time  step  to  the  particles.  Particles  are  given  a  specified  tetrahedral 
shape  whose  base  area  is  equal  to  four  computational  cells;  however,  there  are  many 
particles  in  each  ceU.  The  height  of  each  particle  is  fixed  and  may  be  constant  for 
all  particles  or  may  vary  from  particle  to  particle.  In  either  case  criteria  are  estab¬ 
lished  for  the  number  of  particles  required  for  each  layer.  The  efficacy  of  the  model 
is  illustrated  by  comparing  solutions  with  those  from  an  exact  solution  of  a  nonlinear 
reduced  gravity  model  of  a  parabolic  lens.  The  particle-in-cell  model  reproduces  the 
essential  characteristics  of  the  reduced  gravity  model  including  exceptional  resolution 
of  the  time  varying  surface  front  of  the  lens. 
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1  INTRODUCTION 

Most  simulations  of  classical  fields  (electricity  and  magnetism,  elasticity,  fluid  dynamics) 
have  utilized  Eulerian  or  fixed  grid/element  techniques.  But  as  noted  by  Olim  [1],  time  steps 
utilized  by  these  schemes  may  be  restricted  by  stability  or  damping  so  as  to  be  considerably 
smaller  than  the  time  interval  appropriate  for  the  system  evolutionary  time  scale.  As  a 
consequence,  Lagrangian  or  particle  techniques  have  become  fashionable  for  problems  with 
laxge  velocities. 

With  pure  Lagrangian  techniques,  the  particles  comprise  the  grid.  These  methods  have 
long  been  popular  with  vortex  models  since  the  early  work  of  Christiansen  [2]  and  Chris¬ 
tiansen  and  Zabusky  [3].  Some  recent  examples  are  Dukowicz  and  Meltz  [4],  Winckelmans 
and  Leonard  [5]  and  Russo  [6].  However,  as  pointed  out  in  [1],  a  potential  problem  with  these 
methods  is  the  deterioration  of  spatial  resolution  in  those  parts  of  the  calculation  domain 
with  low  vortex  density.  Also,  with  this  method  it  is  necessary  to  calculate  a  nontrivial 
interaction  of  each  particle  with  all  other  particles  at  each  time  step  using  the  Biot-Savart 
law.  Such  calculations  increase  greater  than  linearly  with  the  number  of  particles  and  are 
particularly  susceptible  to  round-off  errors. 

Semi-Lagrangian  or  particle-in-cell  (PIC)  schemes  can  overcome  some  of  these  prob¬ 
lems.  With  this  method  all  gradient-type  terms  in  the  conservation  of  mass  and  momentum 
equations  are  computed  at  fixed  grid  points  while  the  material  derivatives  are  computed  at 
particles.  Much  of  the  calculation  is  devoted  to  interpolating  particle  properties  to  the  fixed 
grid,  calculating  gradients,  and  then  interpolating  gradient  values  from  the  grid  back  to  the 
particles.  As  with  vortex  methods  there  is  a  trade-off  between  resolution  and  number  of 
particles.  The  computational  load  for  large  particle  numbers,  however,  does  not  appear  to 
be  as  severe  with  PIC  schemes  as  it  does  with  vortex  models  since  it  increases  only  linearly 
with  the  number  of  particles.  This  method  was  first  developed  by  Harlow  [7]  for  fluids.  Since 
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then  there  have  been  extensive  applications  by  Brunei  et  al.  [8],  Brackbill  and  Ruppel  [9] 
and  O’Rourke  et  al.  [10]  to  name  but  a  few. 

The  concern  here  is  with  geophysical  fluid  dynamics.  This  discipline  poses  two  difficult 
problems  for  modelers.  One  is  the  enormous  range  of  scales  and  consequent  phenomenology 
for  which  the  model  must  account.  With  basin  scale  primitive  equation  Eulerian  models  the 
grid  spacing  generally  is  much  too  coarse  to  resolve  all  the  energetic  motions.  The  other 
problem  is  fronts  or  interfaces.  The  formation  and  evolution  of  many  types  of  atmospheric 
or  oceanographic  fronts  have  no  analog  in  fluid  mechanics  other  than  multi-fluid  systems.  As 
noted  in  [10],  PIC  methods  which  allow  for  conservation  of  a  material  property  at  particles 
can  be  very  effective  for  tracking  material  interfaces.  Despite  their  advantages  vis-a-vis  Eu¬ 
lerian  methods  there  have  been  very  few  applications  of  Lagrangian  methods  in  geophysical 
fluid  dynamics.  Notable  exceptions  are  Zabusky  and  McWilliams  [11]  and  Hooker  [12],  who 
used  point-vortex  models  to  study  geostrophic  vortices,  and  Pavia  and  Cushman-Roisin  [13, 
14],  Pavia  [15]  and  Mathias  [16],  who  used  PIC  methods  to  study  ocean  fronts  and  merging 

of  ocean  eddies. 

The  PIC  applications  in  geophysical  fluid  dynamics  have  been  used  in  layer  models.  The 
model  equations  result  from  a  vertical  integration  of  the  hydrodynamic  equations  along  with 
an  appeal  to  quasi-hydrostatic  equilibrium.  If  the  height  of  individual  particles  is  specified 
either  a  priori  or  by  some  pseudo-equation  of  state  then  there  is  no  need  to  calculate  the 
horizontal  divergence  and  solve  the  conservation  of  mass  equation  at  each  particle.  This  is 
considerably  simpler  than  the  standard  Eulerian  approach  of  solving  a  Poisson  problem  for 
the  pressure  at  each  time  step.  However,  there  is  still  the  problem  of  interpolating  particle 
heights  to  the  grid  and  the  gradients  back  to  the  particles  in  a  self-consistent  manner.  A 
second  issue  is  the  number  of  particles  to  use.  To  date,  this  has  been  empirically  set  as  an 
ad  hoc  balance  between  accuracy  and  computational  feasibility. 

There  are  several  unusual  features  about  our  approach.  First,  we  use  the  full  shallow 
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water  hydrodynamic  equations.  Second,  we  use  an  interpolation  routine  that  exactly  kills 
particle  self-motion  caused  by  fictitious  pressure  gradients  associated  with  the  distribution 
of  the  particle  heights  on  the  grid.  Moreover,  our  approach  does  not  require  smoothing  or 
damping  other  than  that  resulting  from  the  interpolation  routines.  Finally,  we  are  able  to 
relate  the  number  of  particles  to  a  desired  spatial  resolution.  This  is  achieved  by  requiring 
the  total  volume  of  particles  to  be  the  same  as  the  fluid  volume.  Since  the  particle  volume 
is  fixed  our  method  guarantees  conservation  of  volume. 

2  THEORY 

2.1  Two-layer  Model 

The  dynamical  basis  for  the  applications  considered  here  are  the  hydrodynamic  equations  for 
a  two-layer  fluid  in  a  steadily  rotating  coordinate  system.  Following  the  notation  of  Hurlburt 


and  Thompson  [17]  these  are 

-^vi  -f  k  X  /vi  =  — 5'V(/ii  4-  /i2  +  T>)  +  ^  "  +  AVVi,  (1) 

at  Pi^i 

^  +  A.V.v.=0,  (2) 

do 

.  -^V2  -h  k  X  /v2  =  -gy{hi  -f  /i2  +  T>)  +  -)-  ^  ^  +  AVV2,  (3) 

at 

^  +  /l2V-V2  =  0,  (4) 


In  these  equations,  Vj  is  the  horizontal  velocity  vector  for  layer  i;  k  is  the  unit  vertical 
vector;  /  is  Coriolis;  g  is  the  gravity  constant;  hi  are  the  instantaneous  layer  thicknesses 
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relative  to  a  common  datum;  D  is  the  bottom  topography;  tj  are  the  stresses  at  the  surfaces 
j  -  g  (upper  surface),  I  (interface  between  the  layers),  and  B  (bottom  surface);  pi  are  the 
constant  densities  of  each  layer;  A  is  the  horizontal  turbulent  viscosity;  and  g,  =  g{p2-pi)l P2 
is  “reduced”  gravity.  These  equations  are  obtained  by  vertically  integrating  the  instantaneous 
three-dimensional  hydrodynamic  equations.  Figure  1  is  a  cartoon  of  the  geometry  for  the 
problem  considered  below. 

As  the  primary  purpose  here  is  to  test  the  efficacy  of  the  PIC  method  with  special  cases 
it  will  not  be  necessary  to  consider  either  the  surface  stresses  n  or  turbulent  viscosity  A. 
Consequently,  these  terms  will  be  neglected  below.  Furthermore,  since  the  test  problems 
studied  here  are  for  the  “reduced”  gravity  special  case  of  (l)-(4)  it  is  sufficient  to  set  D  =  0. 

The  physical  basis  for  the  reduced  gravity  case  is  that,  to  a  first  approximation,  the 
lower  layer  only  adjusts  hydrostatically  to  accelerations.  Hence,  there  is  negligible  flow  in 
the  lower  layer  and  its  thickness  gradient  V/12  can  be  eliminated  in  favor  of  Vh  from  (3) 
and  substituted  into  (1).  The  result  of  this  manipulation  is 

k  X /v  = (6) 

at 

^  +  /iV  ■  V  =  0.  (7) 

dt 

As  there  is  flow  in  only  the  upper  layer,  the  subscripts  denoting  layers  are  superfluous; 

consequently,  they  were  dropped  in  (6)  and  (7)  and  below. 

At  this  point  it  is  appropriate  to  nondimensionalize  (6)  and  (7).  The  scheme  used  by 
Cushman-Roisin  et  al.  [18]  is  appropriate.  This  scales  hhy  H  (representative  layer  thickness), 
V  by  time  by  f-\  and  v  by  This  scaling  removes  all  explicit  parameter 

dependence.  The  resulting  equations  are 
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^  +  kxv  +  V/i  =  0, 

dt 


(8) 

^  +  hV-v  =  0.  (9) 

dtt 

Since  we  use  a  particle  method  these  last  three  equations  must  be  augmented  by  two  path 
equations: 

^_v  =  0.  (10) 

dt 

Equations  (8)-(10)  constitute  five  nonlinear  coupled  equations  governing  the  flow  in  the 
upper  layer.  We  solve  these  as  an  initial  value  problem  using  the  PIC  technique  described 
below. 

2.2  PIC  Paradigm 

The  basis  of  PIC  methods  is  to  solve  the  d/dt  terms  of  (8)-(10)  at  the  particles  and  the 
gradient  and  divergence  terms  on  a  discrete  grid.  Interpolation  between  the  particles  and 
field  occurs  at  each  time  step.  The  central  assumption  of  the  PIC  method  as  applied  to 
these  geophysical  fluid  dynamics  problems  is  that  the  heights  and  volumes  of  each  individual 
particle  are  invariant  with  position  and  time.  Thus  there  is  no  need  to  solve  (9)  at  individual 
particles.  As  suggested  in  [10],  keeping  each  particle  height  fixed  provides  a  means  for 
accurate  tracking  of  the  interface  between  the  active  and  inert  layers. 

The  height  gradient  appearing  in  (8)  is  found  by  computing  (interpolating)  the  h  field 
to  the  grid  using  the  positions  and  heights  of  the  particles,  computing  V/i  on  the  grid,  and 
finally  computing  (interpolating)  V/t  back  to  each  particle.  In  order  to  clarify  these  ideas 
we  outline  the  algorithm.  Assume  that  the  position,  velocity  and  height  of  each  particle  is 
known  at  time  t.  The  steps  in  the  algorithm  to  advance  these  data  to  time  t  St  are 
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1.  interpolate  the  heights  of  the  particles  to  the  grid, 

2.  calculate  the  finite  difference  approximation  to  V/i  on  the  grid, 

3.  interpolate  V/i  from  the  grid  to  the  particles, 

4.  using  an  appropriate  time  integrator,  simultaneously  integrate  (8)  and  (10)  from  t  to 
t  +  6t.  This  yields  the  position  and  velocity  of  each  particle  at  the  new  time. 

For  geophysical  fluid  dynamics  problems,  the  critical  issue  is  the  calculation  of  the  height 
or  pressure  gradients.  In  the  standard  approach  used  in  [13  -  16],  the  height  field  is  specified 
as  the  sum  of  heights  over  all  particles.  Individual  particle  heights  are  usually  fixed  and  their 

base  areas  are  no  larger  than  a  computation  cell. 

An  alternate  approach  is  employed  here.  There  are  some  unique  aspects  about  our 
method  so  some  explanation  is  offered.  As  in  prior  approaches,  each  particle  is  assigned  a 
finite  volume  where  the  total  volume  of  the  particles  is  the  same  as  the  initial  volume  of  the 
active  layer.  We  differ  from  other  investigations  in  that  we  require  the  base  of  each  particle 
to  be  a  square  with  sides  of  length  2A,  where  A  is  the  fixed  distance  between  cell  centers. 
Thus,  a  typical  particle  base  has  the  area  of  four  cells  and  generally  will  overlap  into  nine 

cells. 

The  height  equation  for  a  particle  is 

Zp  =  (11) 

where  Q  are  the  component  distances  of  the  apex  from  the  nearest  cell  center.  Note  that 
l^^l,  |^j,|  <  A.  The  centerline  height  of  a  particle  is  hp,  which  may  vary  from  particle  to  par¬ 
ticle.  As  can  be  seen  from  (11),  the  height  decays  linearly  to  zero  at  the  particle  boundaries. 
This  shape  was  utilized  by  Hockney  and  Eastwood  [19]  who  termed  it  a  “triangular- shaped 
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cloud.”  Figure  2  is  a  three-dimensional  perspective  of  this  shape.  Within  any  one  grid  cell 
there  will  be  considerable  overlap  of  particle  bases. 

It  must  be  understood  that  the  particles  axe  not  to  be  thought  of  as  solid  particles  that 
cannot  be  interpenetrated  but  as  a  computational  representation  of  fluid  particles.  If  several 
particles  were  to  be  located  at  exactly  the  same  (x,  j/)  location  we  interpret  this  as  meaning 
that  the  height,  h,  at  that  location  is  the  sum  of  the  height  of  the  individual  particles.  These 
particles  do  not  lose  their  identity  as  they  interpenetrate  and  pass  through  each  other  and 
so  the  number  of  particles  is  invariant.  There  are  no  particle-particle  interactions  other  than 

what  arise  through  the  gradient  of  the  height  fleld. 

The  key  to  our  approach  is  the  method  by  which  the  height  of  the  particles  is  apportioned 
to  the  centers  of  cells  which  they  overlie,  the  method  used  to  calculate  height  gradients  on 
the  fixed  grid  composed  of  cell  centers,  and  then  the  interpolation  routine  used  to  move  these 
back  to  the  particles.  Our  experience  has  been  that  these  operations  are  not  independent  but 
must  be  done  in  an  internally  consistent  fashion.  The  principle  used  here  is  to  require  that 
gradients  interpolated  to  a  particular  particle  be  independent  of  the  height  of  that  particle. 
This  insures  no  self-generated  motion. 

Figure  3  illustrates  a  typical  cross  section  along  the  x  axis.  Prom  (11),  the  normalized 
cross-sectional  areas  of  the  particle  in  each  of  the  enveloped  cells  along  this  axis  are 

wj.i  =  (1/2  -  f./A)V2, 

=  {1  -  (l/2)(l/2  +  2(6/A)’]}, 

=  (1/2  +  6/A)V2.  (12) 

where  the  normalization  factor  is  the  cross-sectional  area  of  the  particle.  The  interpolation 
of  the  height  of  this  particle  to  the  overlapped  cell  centers  is  then  given  by  hpWj-i^  hpWj, 
hpWj+i  for  the  ;  -  1,  j,  and  j  -fl  cells,  respectively.  The  hatched,  open,  and  cross-hatched 
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areas  in  figure  3  show  how  the  height  is  partitioned  at  the  three  grid  points. 

The  above  results  apply  to  a  typical  cross  section  for  a  single  particle.  The  two-dimensional 
case  for  a  generic  particle  p  is  achieved  by  prescribing  the  weights  as 

To  summarize,  the  height  at  the  ith,  jth  cell  center  is  the  weighted  average  of  all  particles 
enveloping  the  cell.  This  is  given  by 

Nij 

p=l 

where  p  is  the  particle  number  in  the  cell  and  Nij  is  the  number  of  particles  that  overlap 
into  the  cell  ij.  After  the  heights  are  accumulated  at  all  grid  points,  second  order  accurate 
gradients  are  calculated  and  then  interpolated  back  to  the  particles  using  the  same  weights. 
From  (14),  the  values  of  the  gradients  at  the  nine  relevant  cell  centers  are  readily  expressed 

as 


Di+a,j-i-p  =  [hi+a+l,j+P  ~ 
for  the  X  component  of  the  gradient  and 

=  [hi+a,j+P+l  ~  ^i+a,  P  ~ 

for  the  y  component. 

Interpolation  of  these  gradients  back  to  a  particular  particle  P  in  cell  i,  j  is  independent 
of  hp.  To  see  this,  focus  on  the  x  component,  (15A).  The  interpolation  of  this  component 

back  to  P  is 
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(16) 


1 

=  E  U^P(t+a,;+l)A+a,  ;+l  + 

or=  — 1 

1  1 

^P(i+aJ)^i+a,j  -i-  ^i^P(t+Of,j-l)A+Qr,i-l- 

a=^l  a=-l 

Using  (14)  and  (15A)  the  first  sum  in  (16)  is 


2A  '^p(i+oi,j}-C^i+a,j+l 

Of=  — 1 

^P(i-I,i+1)[  £  ^php-  ^  +  X)  Wphp-  Wphp] 

p=i  p=i  p=i  p=^ 

•/Vi+2,i+i 

+i^p(t+i„i+i)[  X)  '^p^p  ~  X^  i"p^p]- 

p=l  p=l 

Now,  hp  does  not  appear  in  the  Ni-2yj+i  ^i+2j+i  sums  since  P  is  not  in  either  cell.  It 
also  appears  once  in  each  of  the  other  four  sums;  however,  the  terms  exactly  cancel.  It  is 
readily  seen  that  this  is  true  for  each  of  the  other  sums  in  (16)  as  well  as  the  y  component 
of  the  gradient. 

The  cancelation  of  particle  height  in  the  interpolated  gradient  to  a  particle  is  analogous 
to  the  condition  that  second  order  accurate  finite  differences  at  grid  points  do  not  depend 
on  the  value  of  the  function  at  those  grid  points.  The  importance  of  this  to  the  present 
application  is  that  there  is  no  self-induced  motion  of  a  particle.  This  result  arises  for  two 
reasons.  First,  the  same  weights  are  used  for  the  interpolation  from  particles  to  the  grid 
and  then  back  to  the  particles  and  not  because  of  a  specific  form  such  as  (12).  Second,  the 
order  of  accuracy  of  the  derivative  (second  order  in  our  case)  is  not  higher  than  that  of  the 
interpolation  routine.  Fourth  order  accurate  derivatives  used  in  conjunction  with  a  lower 
order  interpolation  scheme  may  produce  self  motion  and  thus  would  not  be  as  accurate  as 
the  second  order  schemes. 
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This  method  has  several  other  attributes.  First,  the  weights  are  computationally  efScient 
since  they  are  not  difficult  to  calculate  and  axe  used  both  for  interpolation  from  particles  to 
grid  points  and  back  to  the  particles.  Second,  (9)  is  automatically  satisfied  since  the  heights 
and  volumes  of  the  particles  are  assigned  initially  and  axe  fixed.  Moreover,  it  is  not  necessary 
to  interpolate  particle  velocities  as  long  as  there  is  no  viscosity  in  the  problem.  Finally,  any 
stable  time  integrator  can  be  used  to  advance  the  position  each  particle. 

The  major  tradeoff  in  choosing  the  integrator  is  that  between  accuracy  and  computational 
work.  For  example,  a  two  step  integrator,  such  as  a  second  order  Runge-Kutta  scheme,  has  an 
accuracy  of  0{(StY)  but  requires  two  evaluations  of  V/i  at  each  particle.  Thus  this  integrator 
would  require  two  interpolations  of  the  particle  heights  to  the  grid,  two  applications  of  the 
finite  difference  form  of  V  to  ^  on  the  grid  and  two  interpolations  of  V/i  to  the  particles  for 
each  position  advancement.  On  the  other  hand  a  first  order  accurate  integrator  with  a  much 
smaller  time  step  could  provide  the  same  degree  of  accuracy  with  less  computation  per  time 
step. 

2.3  Initialization  of  Height  Field 

We  shall  illustrate  this  method  in  the  next  section  by  an  application  to  an  initially  circular 
lens.  For  initial  conditions  we  choose  a  parabolic  lens  that  has  been  the  subject  of  consid¬ 
erable  study  [18,  20  -  35].  In  these  studies  the  lens  velocity  and  thickness  were  specified 

as 


h  =  ho(t)  -f  Bij(t)xiXj,  for  h>0 
Vi  =  Gij(t)xj,  for  A  >  0 

Vi  =  0,  elsewhere  (17) 

Here  the  summation  convention  is  used  for  repeated  indices  and  the  spatial  coordinate  Xi  is 
measured  from  the  center  of  mass  of  the  lens.  Note  that  the  velocity  field  is  discontinuous 
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at  /i  =  0  as  is  V/i.  Moreover,  the  discontinuous  frontal  boundary,  h  =  0  must  be  calculated 
as  part  of  the  solution.  Conventional  gridded  methods  are  not  well  suited  for  problems  of 
this  sort.  " 

Substitution  of  (17)  into  (8)  and  (9)  yields  8  coupled  nonlinear  ordinary  differential 
equations  for  ho{t)  and  %(i).  As  shown  in  [24]  and  [34]  the  resulting  solutions 

can  be  represented  as  a  nonlinear  superposition  of  rotational,  deformational  and  horizontal 
divergence  modes.  The  lens  frontal  boundary  generally  may  oscillate  between  a  circle  and  an 
ellipse  which  tends  to  rotate  anticyclonically.  Particle  motion  can  be  quite  complicated;  see 
for  example  sample  trajectories  in  [31].  Solutions  based  on  (17)  will  be  called  lens  solutions. 
The  nondimensional  equation  for  the  initial  lens  thickness  is 

h  =  ho[l  -  (r/R)%  (18) 

where  ho  is  the  centerline  thickness,  r  is  the  radial  coordinate  from  the  lens  center  and 
R  =  [\Bn  +  .^221/2]“^  is  the  radius  of  the  lens.  The  volume  occupied  by  the  lens  is 

Vl  =  27r  r  rhdr  =  (7r/2)/i,i7".  (19) 

Jo 

Integration  of  (11)  shows  the  volume  of  a  particular  particle  to  be  Vp  =  A'^hp  so  the  total 
volume  occupied  by  the  particles  is 

l/r  =  A^E'*p,  PO) 

P  =  1 

where  N  is  the  total  number  of  particles.  Clearly,  Vt  =  Vl  must  be  required. 

Consider  first  the  case  where  hp  is  the  same  for  all  particles.  Then, 

N  =  Vt/Vp  =  {i:l2){holhp){RlAf.  (21) 

To  estimate  N,  divide  the  lens  into  concentric  circular  annuli  surrounding  an  axial  cylin¬ 
der  centered  at  the  lens  center.  The  radial  width  of  an  annulus  is  2A/ d  with  d  >  2,  and  the 
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height  of  the  lens  along  the  axis  of  an  annulus  is  given  by  (18).  In  the  limiting  case,  d  =  2, 
the  annulus  width  is  the  same  as  the  computational  cell  discussed  in  the  previous  section. 
Since  the  amnuli  are  used  for  the  initial  distribution  of  particles  in  the  lens  it  is  appropriate 
to  have  the  annuli  spacing  finer  than  that  of  the  computational  cells. 

The  axis  of  the  outermost  annulus  is  at  radius  r  =  i?  -  A/d.  Then  the  height  along  this 

axis  is 


hL  =  ho{l^lRd){2-^IRd).  (22) 

Obviously,  hp  <  with  the  equality  holding  for  just  one  particle  in  each  cell  on  the  lens 
boundary.  Using  this  limiting  case  for  hp  in  (21)  gives 

N  =  (7r/2)(E/A)"d/(2  -  A/Rd).  (23) 

All  having  the  same  height,  these  particles  must  be  distributed  nonuniformly  in  the  lens. 
We  shall  use  R  =  10"S  A  =  4  x  10"^  and  d  =  16  so  (23)  suggests  N  ~  0(2.5  x  10®)  for 
/jp  ^  However,  requiring  ten  particles  in  each  cell  in  the  outer  boundary  increases  N  by 

an  order  of  magnitude. 

The  number  of  particles  can  be  reduced  by  decreasing  the  horizontal  resolution,  i.e., 
increasing  A/R.  From  (23)  it  is  seen  that  changes  in  resolution  can  cause  dramatic  changes 
in  the  particle  count.  An  order  of  magnitude  change  in  resolution  produces  three  orders  of 

magnitude  change  in  the  number  of  particles. 

One  of  the  problems  that  arises  with  a  constant  hp  is  that  the  number  of  particles  in  each 
cell  will  vary  with  radius  with  fewer  particles  in  the  outermost  cells.  This  can  be  overcome 
by  allowing  hp  to  vary  with  radius  as  well.  The  volume  of  an  annulus  centered  at  r  is 

VA  =  2ir  f  ^  ^  rhdr  =  ATrhor{A/d)[l  -  {r/R)^  -  {A/Rd)^]. 

Jr-A/d 
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The  number  of  particles  required  to  fill  this  volume  is 

Na  =  Va/Vp  =  ATT[ho/hp{r)]{r/Ad)[l  -  {r/Rf  -  {A/ Rdf],  (24) 

The  number  of  cells  in  this  annulus  is  the  ratio  of  the  annulus  area  to  cell  area.  A  simple 
calculation  gives  this  number  as  Nc  =  A-KrlAd.  If  the  number  of  particles  in  each  cell, 
N:jj,  =  N a! Nc,  is  to  be  constant  the  particles  heights  must  be  distributed  as 

hp{r)  =  (/io/iV#)[l  -  {rlRf  -  {NRdf]-  (25) 

Adjusting  the  height  of  the  particles  to  keep  the  number  of  particles  in  a  cell  constant 
also  changes  the  requirements  on  the  total  number  of  particles.  The  cross-sectional  area 
of  the  lens  annulus  is  irlR^  —  (A/d)^].  Thus,  the  number  of  cells  needed  to  cover  the  lens 
annulus  is 

Nc.n  =  ir[R^  -  (A/d)2]/A^  (26) 

and  tlie  total  number  of  particles  in  the  lens  annulus  is 

Nta  =  IV#  •  -A^ceu.  (27). 

To  complete  this  analysis  it  is  necessary  to  fill  the  inner  cylindrical  core.  Its  volume  is 

Vi  =  2Trf  r^dr  =  7r/i<,(A/d)^[l  —  (1/2)(A/J?d)^].  (28) 

Jo 

The  number  of  particles  required  to  fill  this  void  is 

N,  =  -  (l/2)(A/iM)*]. 

This  is  determined  by  an  arbitrary  choice  of  hp(/).  Here  we  use  Ni  ~  10®. 


(29) 


A  critical  issue  with  this  approach  is  an  appropriate  number  of  particles  for  each  cell 
N^.  Using  different  particle  geometry  than  we  do,  other  investigators  [13  -  15]  used  about 
15  particles  per  cell.  It  should  be  noted  that  these  investigators  used  smoothing  at  every 
timestep.  With  our  geometry  and  elimination  of  smoothing  it  appears  that  ~  5  X  10  is 
required  for  reasonable  resolution  and  accuracy. 

With  5  •  10^  the  total  number  of  particles  is 

Nt  =  Nta  +  Ni  10^.  (30) 

Finally,  it  should  be  noted  that  we  use  a  polar  coordinate  system  here  only  to  initialize 
the  distribution  of  particles.  The  calculations  presented  in  the  next  section  are  done  with  a 
rectangular  grid  even  for  those  problems  which  exhibit  radial  symmetry. 

3  APPLICATION  TO  AN  ISOLATED  LENS 

The  first  simulation  discussed  here  was  designed  to  test  the  ability  of  the  PIC  ansatz  to 
reproduce,  with  small  phase  error  and  amplitude  distortion,  a  single  frequency  in  a  nonlinear 
fiow  and  to  provide  high  resolution  of  the  associated  oscillating  front.  Here  and  in  the  next 
problem  we  shall  use  the  variable  height  formulation.  It  is  stressed  that  there  is  no  smoothing 
in  the  results  presented  below  other  than  what  occurs  with  the  interpolation. 

As  a  control  for  the  PIC  simulation,  comparisons  will  be  made  with  numerical  solutions 
of  the  lens  model.  In  the  latter  model  the  height  field  is  parabolic  and  the  velocity  field 
linear  with  respect  to  the  distance  from  the  lens  center  for  all  time.  These  constraints  are 
imposed  only  initially  on  the  PIC  model.  In  both  cases  a  second  order  accurate  Runge-Kutta 
integrator  was  used. 

The  initial  velocity  field  for  the  PIC  model  was  taken  as  the  same  as  the  lens  model,  i.e., 
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ui(0)  =  (G/2)xi  -  {Gr)x2, 

^;2(0)  =  (C?fi)a:i  -  {G/2)x2,  (31) 

where  G  =  Gu{0)  +  (?22(0)  is  the  initial  horizontal  divergence  and  Gr  =  [C?i2(0)  -  (?22(0)]/2 
is  the  initial  spin.  Numerical  values  for  this  experiment  are  given  in  the  table. 

For  the  lens  model  [21]  obtained  an  analytic  solution  for  ho,  Bij  and  Gij,  which  is  called 
here  the  pulson  solution.  The  nondiagonal  components  of  Bij  and  the  symmetric  component 
of  Gij  are  zero  for  all  time  while  the  remaining  components  as  well  as  ho  oscillate  at  the 
Coriolis  frequency  f .  Unlike  most  nonlinear  problems  the  oscillation  frequency  of  the  pulson 
is  independent  of  the  initial  conditions.  Of  particular  interest  here  are  the  analytic  expres¬ 
sions  for  the  centerline  height,  ho,  and  [|J5ii  4-  JB22I/2]  ^  since  these  two  variables  determine 
the  height  field  and  the  evolution  of  the  lens  boundary.  These  are  found  to  be 

K{t)  = 

h{t)  =  K{t)  -  (AB/4)r"(i),  (32) 

where  if  is  the  layer  thickness  scale  at  (  =  0,  Ab  is  given  by  the  initial  conditions  of  Bn + S22 
and 


Tit)  =  [A +  'f sin ft]-\  A>\^\. 

Also,  A  and  7  axe  given  by  initial  conditions  of  Gij. 

Figure  4  shows  typical  trajectories  for  the  PIC  (solid  curve)  and  lens  solutions  (dashed 
curve).  Both  paxticles  start  at  same  initial  position  shown  as  a  triangle  and  both  exhibit  the 
same  cycloidal  behavior  characteristic  of  this  solution.  The  boxes  show  the  position  at  one 
day  intervals.  This  motion  is  composed  of  general  anticyclonic  rotation  about  the  lens  center 
with  small  loops  at  the  inertial  frequency.  As  seen  in  panel  a,  the  trajectories  are  virtually 
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identical  for  the  first  1.5  days;  however,  by  day  2  some  offset  is  apparent.  Although  there  is 
some  growth  in  the  offset  it  is  mostly  completed  by  day  5  as  seen  in  the  starting  and  ending 
offsets  in  panel  b.  The  PIC  solution  exhibits  some  irregulaxity  in  the  spacing  of  the  inertial 
loops.  The  period,  however,  is  very  nearly  inertial  as  expected. 

Figure  5  compares  the  PIC  solution  for  the  centerline  height  ho  (solid  curves)  with  the 
analytic  solution  given  by  (32).  This  figure  shows  that  the  PIC  solution  accurately  reproduces 
both  the  amplitude  and  phase  of  the  analytic  solution.  Figure  6  compaxes  the  PIC  height 
field  after  10  days  (solid  curve)  with  the  initial  height  field.  Because  the  solution  is  periodic, 
these  fields  should  be  identical.  Again,  the  comparison  is  excellent  although  there  is  some 
indication  that  near  the  lens  center  the  PIC  simulation  is  slightly  thinner  than  the  exact 
solution.  This  is  also  reflected  in  the  slight  increase  in  the  simulation  thickness  at  mid  radii. 
The  difference  between  the  initialized  field  for  PIC  and  the  analytic  solution  could  not  be 
seen  at  normal  plotting  scales. 

The  above  flow  field  had  vorticity  and  horizontal  divergence  but  no  deformation.  To  test 
the  efficacy  of  the  PIC  model  in  the  presence  of  a  time  dependent  deformation  component 
of  the  velocity  field,  the  initial  velocity  (31)  was  replaced  by 

ui(0)  =  (Gj2  +  Gj^)xi  +  (t?jv  Gii)x2, 

1)2(0)  =  (Gn  +  Gr)xi  +  {G/2  —  Gn)x2-  (^^) 

Here  Gn  =  [(?ii(0)-(?22(0)]/2  is  the  initial  normal  deformation  and  Gs  =  [C?i2(0)+C?2i(0)]/2 
is  the  initial  shear  deformation.  Numerical  values  axe  given  in  the  table. 

The  solution  for  the  lens  case  was  obtained  by  numerically  integrating  the  eight  equations 
for  ho,  Bij  and  Gij.  See  [31]  and  [33]  for  details  regarding  the  solution  procedure.  This 
solution  shows  that  the  lens  boundary  starts  as  a  circle,  quickly  deforms  into  an  ellipse 


19 


which  rotates  anticyclonically  for  a  brief  period  then  abruptly  deforms  back  to  a  circle.  This 
cycle  is  repeated;  however,  each  new  repetition  of  the  circular  or  ellipse  phase  does  not  have 
the  same  geometric  characteristics  as  the  previous  phase. 

Figures  7  and  8  show  the  evolution  of  two  example  trajectories.  Figure  7  is  typical  of 
trajectories  for  particles  staxted  near  the  lens  center.  This  figure  shows  that  the  trajectory 
for  the  lens  model  particle  tends  to  stay  near  the  center  although  there  is  considerable 
irregularity  in  the  trajectory.  In  contrast  the  PIC  particle  tends  to  wander  away  from  the 
lens  center.  This  is  consistent  with  the  results  of  the  previous  experiments.  Other  than 
exhibiting  anticyclonic  rotation  the  two  trajectories  appear  uncorrelated  at  very  early  time. 

Figure  8  shows  typical  paths  for  particles  started  near  the  initial  lens  boundary.  Here 
the  trajectories  from  the  PIC  and  lens  model  are  in  good  agreement  for  the  first  two  days. 
Thereafter  the  trajectories  exhibit  only  crude  similarities  such  as  anticyclonic  rotation  and 
cycloidal  looping.  As  with  the  previous  experiment,  the  lens  model  trajectory  shows  greater 
excursions  from  the  lens  center  than  does  the  PIC  model. 

Figure  9  compares  the  PIC  centerline  height  (solid  curve)  with  the  numerical  solution  to 
the  lens  model.  As  seen  in  this  figure,  the  solutions  agree  quite  well  through  day  2.  After 
this  time  the  solutions  agree  well  in  phase  but  the  amplitude  of  the  PIC  solution  tends  to 
be  slightly  smaller  than  the  lens  solution.  This  is  consistent  with  the  indication  from  Fig.  7 
that  with  the  PIC  solution  the  lens  tends  to  flatten. 

Figures  10  and  11  compare  the  evolution  of  the  lens  boundary  determined  by  the  PIC 
model  with  the  lens  solution.  The  detailed  evolution  for  the  first  day  is  shown  in  figure 
10  while  figure  11  shows  snapshots  of  the  boundary  over  a  10  day  period.  It  is  seen  in 
these  figures  that  up  to  day  2.5  the  PIC  and  lens  solutions  are  in  close  agreement.  Both 
quickly  evolve  to  an  elongated  ellipse  which  rotates  anticyclonically.  At  day  2.5  the  boundary 
has  evolved  to  a  more  circular  shape.  Thereafter  the  PIC  solution  continues  to  display  an 
elliptical  mode;  however  it  is  not  nearly  as  pronounced  as  that  of  the  lens  solution.  By  day 
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8  there  is  some  indication  of  a  Kelvin-Helmholtz  type  instability  at  the  boundary  of  the  PIC 

solution;  however,  the  growth  rate  does  not  seem  large. 

The  discrepancies  between  the  two  model  simulations  begin  around  day  2.  There  is 
no  indication  of  numerical  instability  in  either  case  and  no  indication  of  a  hydrodynamic 
instability  in  the  PIC  simulations  other  than  a  possible  weak  Kelvin-Helmholtz  instability 
late  in  the  simulation.  Note  that  hydrodynamic  instabilities  are  precluded  by  the  formulation 
of  the  lens  model  since  the  height  and  velocity  profiles  axe  required  for  all  time  to  be  parabolic 
and  linear  respectively.  Thus  either  solution  could  be  representative  of  real  oceanic  lenses. 
Since  the  PIC  model  has  fewer  restrictions  we  prefer  it. 

The  PIC  calculations  described  above  used  approximately  10®  particles.  In  order  to 
illustrate  the  deterioration  in  the  solution  when  less  particles  are  used  the  deformation  case 
was  repeated  with  only  about  10®  particles.  Figure  12  shows  the  evolution  of  the  centerline 
height  for  the  two  cases.  It  is  seen  that  the  peaks  and  troughs  become  ragged  for  the  reduced 
particle  case.  Note  that  for  plotting  purposes  the  10®  particle  case  is  shown  as  the  dotted 
curve  while  the  10®  particle  case  is  shown  as  the  solid  curve. 

4  CONCLUSIONS 

The  results  presented  in  the  previous  section  are  very  encouraging.  The  comparison  of  the 
simulations  and  the  pulson  analytic  solution  showed  that  the  PIC  method  used  here  was  able 
to  reproduce  both  the  amplitude  and  phase  of  the  analytic  pulson  solution  with  negligible 
distortion.  The  case  of  the  deforming  vortex  provided  a  good  test  of  the  method  for  handling 
typical  evolutions  of  fronts  in  which  local  horizontal  divergence,  vorticity  and  deformation  are 
all  important.  The  time  and  space  scales  of  the  front  oscillations  would  be  computationally 
costly  to  treat  by  conventional  gridded  primitive  equation  models.  The  agreement  between 
the  PIC  simulations  and  the  numerical  results  from  the  lens  model  indicate  that  the  former 
is  capable  of  treating  problems  with  dynamically  evolving  interfaces. 
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The  formulation  used  here  allows  for  the  particles  to  interpenetrate  so  that  a  given  spatial 
position  generally  will  be  contained  within  the  volumes  of  many  particles.  This  means  that 
with  this  formulation  it  is  not  appropriate  to  consider  material  properties  as  characteristic 
of  a  single  particle  occupying  a  particular  position  at  a  particular  time.  Rather  material 
properties  should  be  viewed  as  a  suitable  weighted  average  of  all  the  particles  influencing 
the  position  at  the  particular  time. 

It  is  noteworthy  that  the  results  achieved  here  did  not  require  any  smoothing  other  than 
that  implicitly  contained  in  the  interpolation  and  differentiation  operations.  Perhaps  this  is 
due  to  our  choice  of  the  particle  shape  and  weighting.  This  approach  may  be  better  suited 
to  simulating  a  continuous  fluid  media  than  schemes  which  use  simpler  shapes  and  weighting 
functions. 

There  are  still  unresolved  issues  regarding  the  utility  of  this  approach,  however.  No 
one  has  yet  introduced  a  fully  dynamically  interacting  lower  layer.  By  prescribing  a  quasi- 
geostrophic  response  in  the  lower  layer  [16]  was  able  to  use  conventional  gridded  methods  for 
obtaining  the  solution  in  that  region  while  using  PIC  methods  in  the  upper  layer.  Using  PIC 
methods  in  both  layers  should  be  quite  challenging.  Another  question  is  to  determine  how 
well  the  method  treats  streaming  flows.  In  this  case  the  particles  are  advected  out  of  one 
part  of  the  domain  and  must  be  replaced  by  inflow  from  another  part.  There  is  no  experience 
with  PIC  methods  in  a  geophysical  fluid  dynamics  setting  with  problems  of  this  sort.  Finally, 
it  will  be  important  to  consider  surface  stress-driven  flow  and  flow  over  topography.  In  these 
cases  it  will  be  necessary  to  incorporate  viscosity,  which  will  require  some  modifications  in 
the  approach  used  here. 

Acknowledgements.  This  study  was  supported  in  part  by  the  Office  of  Naval  Research 
under  contract  N00014-91-J-1560  and  ONR-AASERT  contract  N00014-93-1-0842.  A.D.  Kir- 
wan,  Jr.  and  C.E.  Grosch  acknowledge  the  Samuel  L.  and  Fay  M.  Slover  endowment  to  Old 


22 


Dominion  University.  The  particle  geometry  used  here  was  suggested  to  us  by  R.  Hockney. 
We  benefited  greatly  from  discussions  with  him  in  September  1994.  Finally,  the  authors 
gratefully  acknowledge  the  technical  assistance  of  K.  Gregory. 


23 


References 


[1]  M.  Olim.  J.  Comput.  Phys.,  112,  253,  1994. 

[2]  J.  P.  Christiajisen.  J.  Comput.  Physics,  13,  363,  1973. 

[3]  J.  P.  Christiansen  and  N.J.  Zabusky.  J.  Fluid  Mech.,  61,  219,  1973. 

[4]  J.  K.  Dukowicz  and  B.  Meltz.  J.  Comput.  Phys.,  99,  115,  1992. 

[5]  G.  S.  Winckelmans  and  A.  Leonard.  J.  Comput.  Phys.,  109,  247,  1993. 

[6]  G.  Russo.  J.  Comput.  Phys.,  108,  84,  1993. 

[7]  F.  H.  Harlow.  Meth.  Comput.  Physics,  3,  319,  1964. 

[8]  F.  Brunei,  J.  N.  Leboeuf,  T.  Tajima,  and  J.  M.  Dawson.  J.  Comput.  Phys.,  43,  28, 
1981. 

[9]  J.  U.  Brackbill  and  H.  M.  Ruppel.  J.  Comput.  Phys.,  65,  314,  1986. 

[10]  P.  J.  O’Rourke,  J.  U.  Brackbill,  and  B.  Larrouturou.  J.  Comput.  Phys.,  109,  37,  1993. 

[11]  N.  J.  Zabusky  and  J.  C.  McWilliams.  Phys.  Fluids,  25,  2175,  1982. 

[12]  S.  B.  Hooker.  PhD  thesis.  University  of  Miami,  Miami,  Florida,  1987. 

[13]  E.  G.  Pavia  and  B.  Cushman- Roisin.  J.  Geophy.  Res.,  93,  3554,  1988. 

[14]  E.  G.  Pavia  and  B.  Cushman- Roisin.  J.  Phys.  Oceanog.,  20,  1886,  1990. 

[15]  E.  G.  Pavia.  PhD  thesis,  Florida  State  University,  Gainesville,  Florida,  December  1989. 

[16]  B.  J.  Mathias.  Master’s  thesis,  Thayer  School  of  Engineering,  Dartmouth  College,  New 
Hampshire,  April  1992. 


24 


[17]  H.  E.  Hurlburt  and  J.  D.  Thompson.  J.  Phys.  Oceanog.,  3,  16,  1973. 

[18]  B.  Cushman-Roisin,  W.  H.  Heil,  and  D.  Nof.  J.  Geophy.  Res.,  90(C6),  11,756,  1985. 

[19]  R.  W.  Hockney  and  J.  W.  Eastwood.  Computer  Simulation  Using  Particles.  Institute 
of  Physics,  1992. 

[20]  G.  R.  Goldsbrough.  In  Proc.  R.  Soc.  London,  A30,  157,  1930. 

[21]  F.  K.  Ball.  J.  Fluid  Mech.,  19,  240,  1963. 

[22]  F.  K.  Ball.  J.  Fluid  Mech.,  22(3),  529,  1965. 

[23]  W.  C.  Thacker.  J.  Fluid  Mech.,  107,  499,  1981. 

[24]  W.  R.  Young.  J.  Fluid  Mech.,  171,  101,  1986. 

[25]  B.  Cushman-Roisin.  Tellus,  39(A),  235,  1987. 

[26]  P.  Ripa.  J.  Fluid  Mech.,  183,  343,  1987. 

[27]  B.  R.  Ruddick.  J.  Phys.  Oceanog.,  17,  741,  1987. 

[28]  C.  Rogers.  Phys.  Lett.  A,  138,  267,  1989. 

[29]  D.  Brickman  and  B.  Ruddick.  J.  Geophy.  Res.,  95,  9657,  1990. 

[30]  A.  D.  Kirwan,  Jr.,  A.  W.  Indest,  J.  Liu,  and  N.  Clark.  J.  Geophy.  Res.,  95,  18,057, 
1990. 

[31]  A.  D.  Kirwan,  Jr.  and  J.  Liu.  In  A.  R.  Osborne,  editor.  Nonlinear  Topics  in  Ocean 
Physics,  Course  CIX,  pages  99-132,  North  Holland,  Amsterdam,  1991.  International 
School  of  Physics,  ‘Enrico  Fermi’. 

[32]  T.  Kawano  and  S.  Nakamoto.  J.  Geophy.  Res.,  97,  12,659,  1992. 


25 


[33]  A.  D.  Kirwan,  Jr.,  B.  L.  Lipphardt,  Jr.,  and  J.  Liu.  Inter.  J.  Eng.  Sci.,  30(10),  1361, 
1992. 

[34]  A.  D.  Kirwan,  Jr.,  B.  L.  Lipphaxdt,  Jr.,  and  K.  L.  Gregory.  In  S.K.  Majumdar,  E.  W. 
Miller,  G.  S.  Forbes,  R.  F.  Scbmalz,  and  A.  A.  Panah,  eds.,  The  Oceans:  Physical- 
Chemical  Dynamics  and  Human  Impact.  The  Pennsylvania  Academy  of  Sciences,  1994. 

[35]  A.  D.  Kirwan,  Jr.  and  B.  L.  Lipphardt,  Jr.  J.  Mar.  Sys.,  4(2-3),  95,  1993. 


26 


Table:  Initial  Conditions 


Pulson  Case 

Deformation  Case 

h/Q 

4.875  X  10-^ 

4.875  X  10-^ 

{Bn 

+  •S22)/2 

-9.75  X  10-^ 

-9.75  X  10-2 

Bi2i 

B21 

0 

0 

Gr 

-0.25 

-0.25 

Gn 

0 

0.2 

Gs 

0 

0.1 

G 

0.6 

0.6 
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Figure  Captions 


Figure  1  Cartoon  of  the  two  layer  model. 

Figure  2  Perspective  picture  of  the  shape  of  a  particle. 

Figure  3  Cross  section  of  a  particle  along  the  x  axis.  The  hatched,  open,  and  cross-hatched 
areas  divided  by  A  give  the  weights  specified  in  (12).  The  two  vertical  dashed  lines 
are  the  cell  boundaries. 

Figure  4  Typical  trajectories  for  the  PIC  (solid  curve)  and  lens  model  (dashed  curve)  for  the 
pulson  case.  The  triangle  is  the  starting  point  for  both  trajectories  and  the  boxes  give 
the  positions  at  one  day  intervals.  Panel  a  shows  the  first  five  days  and  panel  b  the  last 
five  days  of  the  simulation.  The  dotted  line  in  panel  a  gives  the  initial  lens  boundary. 

Figure  5  Comparison  of  the  PIC  centerline  height  with  the  analytic  solution  for  the  pulson  case. 
The  PIC  result  is  the  solid  curve  and  the  analytic  solution  is  the  dashed  curve. 

Figure  6  Comparison  of  h  at  i  =  10  days  from  the  PIC  model  (solid  contours)  with  the  initial 
h  (dashed  contours)  for  the  pulson  case. 

Figure  7  Typical  trajectories  for  the  PIC  (solid  curve)  and  the  lens  model  (dashed  curve)  for 
particles  starting  near  the  origin  for  the  deformation  case.  Panel  a  is  for  the  first  five 
days  and  panel  b  is  for  the  last  five  days  of  the  simulation.  The  symbol  convention  is 

the  same  as  in  Figure  4. 

Figure  8  Typical  trajectories  for  the  PIC  (solid  curve)  and  the  lens  model  (dashed  curve)  for 
particles  starting  near  the  initial  lens  boundary  for  the  deformation  case.  Panel  a  is  for 
the  first  five  days  and  panel  b  is  for  the  last  five  days  of  the  simulation.  The  symbol 
convention  is  the  same  cis  in  Figure  4. 
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Figure  9  Comparison  of  the  PIC  centerline  height  with  the  numerical  solution  to  the  lens  model 
over  10  days  for  the  deformation  case.  The  PIC  solution  is  the  solid  curve  while  the 
lens  solution  is  the  dashed  curve. 

Figure  10  Comparison  of  the  evolution  of  the  lens  boundary  determined  by  the  PIC  model  with 
that  of  the  lens  model  over  1  day  for  the  deformation  case.  The  PIC  solution  is  the 
solid  curve  while  the  lens  solution  is  the  dashed  curve. 

Figure  11  Comparison  of  the  evolution  of  the  lens  boundary  determined  by  the  PIC  model  with 
that  of  the  lens  model  over  10  days  for  the  deformation  case.  The  PIC  solution  is  the 
solid  curve  while  the  lens  solution  is  the  dashed  curve. 

Figure  12  Comparison  of  the  centerline  height  evolution  for  PIC  simulations  with  approximately 
10®  particles  (dotted  curve)  and  10®  particles  (solid  curve)  for  the  deformation  case. 
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